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Abstract 



A new method which allows one to study multiple coherent reflec- 
tion/transmissions by partially transparent interfaces, (e.g., in multi-layer 
mesoscopic structures or grain boundaries in high-Tc's), in the framework of 
the quasiclassical theory of superconductivity is suggested. It is argued that 
in the presence of interfaces, a straight-line trajectory transforms to a simple 
connected 1-dimensional tree (graph) with knots, i.e. the points where the in- 
terface scattering events occur and pieces of the trajectories are coupled. For 
the 2-component trajectory "wave function" which factorizes the Gor'kov ma- 
trix Green's function, a linear boundary condition on the knot is formulated 
for an arbitrary interface, specular or diffusive (in the many channel model). 
From the new boundary condition, we derive: (i) the excitation scattering am- 
plitude for the multi-channel Andr eev/or dinar y reflection/transmission pro- 
cesses; (ii) the boundary conditions for the Riccati equation; (iii) the transfer 
matrix which couples the trajectory Green's function before and after the 
interface scattering. To show the usage of the method, the cases of a film 
separated from a bulk superconductor by a partially transparent interface, 
and a SIS' sandwich with finite thickness layers, are considered. The electric 
current response to the vector potential (the superfluid density ps) with the 
vr phase difference in S and S' is calculated for the sandwich. It is shown 
that the model is very sensitive to imperfection of the SS' interface: the low 
temperature response being paramagnetic {ps < ) in the ideal system case, 
changes its sign and becomes diamagnetic {ps > 0) when the probability of 
reflection is as low clS cl few percent. 



PACS numbers: 74.20.Fg,74.50.+r„74.80.D,74.80F 
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I. INTRODUCTION 



Many important properties of superconductors are related to surfaces and interfaces, the 
Josephson and proximity effects being well-known examples. In recent years, new rich sur- 
face physics has been found in high-T^ oxides after the identification of the d-symmetry of 
the order parameter. On the theoretical side, studying an interface poses certain problems: 
The method of the quasiclassical Green's functions (for a recent review see ^ )which 
is the main tool in the superconductivity theory, cannot be directly applied here since the 
quasiclassical condition is violated by fast change of the potentials on the atomic distances 
in the vicinity of the interface. As shown by Zaitsev |^, the abrupt changes at a specular 
partially transparent interface can be incorporated into a boundary condition for the qua- 
siclassical Green's functions; the condition is a third order equation for the matrix Green's 
function near the interface. Various forms of the boundary condition have been discussed in 
more recent papers 10-0. New difficulties arise when one attempts to describe the coherent 
reflection/transmission by many interfaces, e.g. in a multi-layer mesoscopic structures or 
grain boundaries network in high-Tc's . In this case, Zaitsev's third order boundary condi- 
tion must be satisfied on each interface, and one encounters the problem of solving a system 
of cubic matrix equations. It is not obvious that a solution to the system of equations exists 
and is unique if it exists. Moreover, some authors [|^^ doubt the very applicabihty of the 
quasiclassical scheme in the many interface geometry: They argue that the quasiclassical 
normalization, which is a vital part of the quasiclassical scheme, is not possible in a double 
layer system with partially refiective interface. 

The purpose of the present paper is to re-examine the theory of the interface in the 
quasiclassical description of superconductivity. A new scheme which allows one to incorpo- 
rate specular as well as diffusive interface (s) into the quasiclassical theory is suggested. To 
make the presentation self-contained, we start with a short introduction to the quasiclassical 
theory of superconductivity. 

As first shown by Bardeen, Cooper, and Schrieffer (BCS) |Ty], the phenomenon of su- 
perconductivity can be understood in the framework of a mean-field type scheme where 
the Cooper correlations are introduced through the pair potential A (generally, a function 
of the momentum p) which is related to electron-electron interaction by a self-consistency 
condition. The mean field A may be introduced directly as a kind of Hartree-type potential, 
or it can be derived in the framework of a more sophisticated Eliashberg theory where the 
pair potential comes as the anomalous self-energy in the Gor'kov equations for the Green's 
function. This truly microscopic approach allows one to perform all the normalizations in 
the spirit of the Landau theory of Fermi liquid and to consider superconductors with a strong 
coupling (see Seren and Rainer and references therein.) 

Whatever the method of derivation, the Gor'kov equation for the matrix Green's function 
gives the basis for studying the BCS-type superconductivity. The quasiclassical theory of 
superconductivity offers an approximate simplified scheme of solving the Gor'kov equation. 
To clarify physics behind the approximations, we analyze first the Bogoliubov - de Gennes 
equation |]T0| that is the effective "Schrodinger equation" corresponding to the Gor'kov 
equation (in the weak coupling limit). 

It is well-known that Cooper's pairing in the superconducting state is conveniently de- 
scribed in the language of the electron-hole coherence. On the mean field level, the ground 
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as well as an excited states of the system are products of single particle states, each of them 
a quantum superposition of electron and hole. The electron, tpe-, and hole, tph-, amplitudes in 
the superposition comprise the 2-component single particle wave function, \l/(r,)f:) = (^^)- 



It obeys the Bogoliubov - de Gennes equation fTQ] 



^^a^UJ"l A* -e(p+fA)-f/JUJ 

where ^ {p) = e(p) —fi, e{p) and /i being the electron band energy and the chemical potential, 
respectively, A is the magnetic vector potential; U (r) is the potential energy. The pairing 
potential A, and, in principle, all other potentials must be found self-consistently. 

For future needs we note that in the vicinity of the the Fermi surface C,{Pf) = 0, the 
electron (hole) with the momentum p pp moves with the Fermi velocity v = +{~)-^^\£,=o- 
The particle energy is close to the Fermi energy Ep ~ vpp, and the de-Broglie wave length 
is Xp is of order of Xp ~ h/pp, pp being a typical momentum on the Fermi surface. 

In the superconductors which are good metals in the normal state, the potentials are 
semiclassical (excluding interfaces and disorder which are discussed later) i.e. they are 
slowly varying functions of the coordinate on the scale of the wave length Xp. Indeed, the 
pair potential A changes at the coherence length ~ hv/A, and one estimates the ratio 
Xp/^o as Xp/^o ~ A/Ep. Also, the validity of a semiclassical treatment of magnetic field B 
requires that Xp « Ib, Ib being the magnetic length, Ib = \J^q/B , $o = hc/2e. Since 
superconductivity exists only at -B < Bc2 ~ ^o/io , the ratio Xp/Ib never exceeds A/Ep. 
Seeing that A ~ Tc, the semiclassical conditions Xp/^o , Xp/Ib <^ I are equivalent to the 
requirement that Tc/Ep << 1. In accordance with the Landau theory of Fermi liquid, this 
condition is always is satisfied if the normal state is metallic. 

Most of the physical effects in metals and superconductors (the Hall and thermoelectric 
effects being notable exceptions) can be described in the simplest approximation where all 
the corrections of order T/Ep ~ Tc/Ep are neglected i.e. in the limit Tc/Ep 0. This is 
the approximation where the quasiclassical theory of superconductivity is valid 



Since ~ limit is equivalent to /i — or large mass m ~ Pf/v oo. 

In this limit of quantum mechanics of noninteracting particles, wave packets do not suffer 
quantum broadening and dynamics becomes completely classical: The particle moves along 
a trajectory, position r{t) and momentum p{t) being well defined. Below we analyse how 
the electrons-hole coherence in the superconducting state changes the situation. 

First, we consider in more detail the classical dynamics of the electron and hole separately. 
The Bogoliubov - de Gennes equation where we put A = for the moment, reads 

=[Up- eA) +U]^c , = - (p + eA) + U] (1.2) 

The two equations transform into each other after the substitution t —t and A ^ — A. 
This means that given a solution ipc{'r,t\{A}) corresponding to the vector potential A, the 
function ipd"^: ~^|{~^}) solves the equation for ip^ in the vector potential A. Therefore, 

7/;,(r,t|{A}) = ^,(r,-t|{-A}), (1.3) 
provided ^c{r, t = 0)=Mr,t = 0) 
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If ^ — s> 0, the centre of electron or hole wave packets moves in the r — p space along the 
trajectory specified by the coordinate Ve^hif) and momentum Pe,h(t) as a function of time t. 
The relation between electron- and hole-trajectories can be expressed in the following way. 

Let re(h)(t\{b}) together with Pe(h)(t\{b}) be the trajectory of the electron (hole) in the 
magnetic field b = rot A. From Eq.( |1.3| ) one can conclude that the corresponding classical 
dynamics of electrons and holes are related to each other in the following way: 

r,{t\{b}) = re(-t|{-6}) , PH{t\{b}) = Pe{-t\{-b} (1.4) 

provided the electron and hole trajectories pass through the same point = Vh = Tq and 

Pe = Ph = P0 att = 0. 

One sees from here that if the magnetic field is absent, 6 = 0, or its influence on the 
classical dynamics is negligible, then 

rh{t)=r,{~t) , ph{t) = p,{-t) , (1.5) 

that is the electron and hole move in opposite directions along the same line (path) in the 
r — p space. However, to the extent the magnetic field influences the orbits, the electron 
and hole paths are different [|13|. (Obviously, the role of the magnetic field may play any 



perturbation violating the time reversal symmetry.) 

Now we are in position to analyse how the electron-hole mixing {i.e. A 7^ in Eq.( |l.l| )) 
changes propagation of the wave packets. Consider a wave packet which is initially purely 
electronic {iph = , t = 0), and assume for the moment that Eq.( |1.5| ) is valid. The electron 
moves classically on a trajectory in the r — p space, and provides a source, A*'0e, in the 



equation for ijjh (see Eq. (|1 . 1| ) ) generating a hole wave. Since A is a slowly varying field, the 
source A*(r)'?/'e(r, t) and ipe{fi t) are peaked at the same point of the r — p space. In other 
words, the hole is created at the point of the current position of the electron and with the 
instantaneous electron momentum. Then, by virtue of Eq. ( |1.5| ), the secondary hole moves 
backwards along the path of the primary electron. In turn, the hole creates new electrons 
which move along the same path etc. It is very important that the multiple processes of the 
electron-hole conversion keep the packet on a line in the r — p space which is nothing but 
the classical trajectory. However, the width of the packet along the trajectory grows linearly 
in time oc vt (at times t > h/A) due to the reverse of the velocity under the electron ^ hole 
conversion processes. 

One sees that, the wave packet in a superconductor experiences broadening even in the 
limit h —>■ 0, and, therefore, a quantum description is unavoidable. Nevertheless, the notion 
of the classical trajectory as a line in the r — p space remains meaningful because the 
quantum broadening occurs only along the line. Ultimately, this important feature is due 
to the time reversal symmetry. It holds to the extent Eq.( |1.5|) is accurate i.e. when one can 
neglect the magnetic Lorentz force in the classical dynamics. 

Note the peculiar role of a magnetic field: the difference in the magnetic bending of 
electron and hole trajectories results in the broadening of the coherent electron-hole wave 
packet in the direction transverse to the classical trajectory. At energies ~ A where the 
electron and hole components have comparable weight, the significance of the Lorentz force 



can be estimated [r^| from the ratio ^ < A/Ep where cOc = \eB/mc\ is the cyclotron 
frequency. Since ^ oc ^, one can consistently neglect the Lorentz since the quasiclassical 
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theory is effectively a tlieory of infinitely heavy particles, m — oo as discussed before. It 
seems that in general case the Lorentz force can be incorporated in a theory of supercon- 
ductivity only by a full quantum approach (see, however, Kopnin's quasiclassical theory of 
the Hall effect |jl5|). Sometimes, the magnetic broadening may turn out to be non-crucial, 
e.g. in a spatially homogeneous case, and then certain simplifications may be possible (see 
e.g. 0). 

A more formal and rigourous analysis of electron- hole coherence on classical trajectories 
can be done using a method first suggested by Andreev The stationary state wave 

function is written as \E'(r,t) = '?/;(r)e^P^'^e"«"^* where ip{r) is a slowly varying function 
(provided \E\ <^ Ep). Plugging \l/(r,t) into the Bogoliubov - de Gennes equation Eq. (|l.l| ), 
and using the approximation 

c \i c 

where the small terms of order (Ai?V)^ are neglected, one gets the Andreev equation. Re- 
arranging terms, the Andreev equation may be written in the following form 

-•V.(f,---_,,,.,t))(4>0 (.6, 

where v is the velocity at the point pp of the Fermi surface, Ps denotes ps = —^A and 
for simplicity U = (as it usually the case because of the efficient screening). The most 
important feature here is that the derivative v • V couples the value of the wave function 
only on straight lines in the direction of the velocity v; the lines are the classical trajectories 



when [7 = [|T8|. In this approximation, the quantum coherence exists only along the 
classical trajectories without any coupling between neighbouring paths. These properties 
are in agreement with the qualitative picture of the wave packet spreading along the classical 
trajectory, discussed previously. One may call the envelope function ip in Eq.(|1.6|) the wave 
function on the classical trajectory. 

After this short review of the quasiclassical approximation, our next step is to include 
the interface into the scheme. In this introductory part of the paper, we present main ideas 
using the language of the wave functions on classical trajectory; a more general approach of 
2-point trajectory Green's function is presented in Section ^ 

The reflection/transmission on an isolated interface (a specular one, to begin with) mixes 
together semi- infinite pieces of classical trajectories (see Fig.p. Each of the pieces is charac- 
terized by the Fermi surface momentum pp, and the corresponding velocity is v; the arrows 
indicate the direction of the velocity. On pieces of trajectories 1 and 2 the velocity is directed 
towards the interface, and we call them in-coming trajectories (or channels); correspond- 
ingly, 1' and 2' are outgoing (pieces of) trajectories. Throughout the paper, the out-going 
"channels", alias for "trajectory", are marked by "prime". 

Note that the in/out classification of the trajectories in accordance with the direction 
of the Fermi surface velocity is unique but it is arbitrary because the electron and hole 
belonging to same channel have the opposite directions of their velocities. For instance, the 
electron coming to the interface on via e.g. the channel 2 (see Fig. |I|) may go away as the 
electron on trajectories 1' and 2' as well as a hole along nominally in-coming trajectory 1. 
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We will call "knot" the region inside of which scattering occurs and the pieces of the 



classical trajectories get "tied" together on the interface ||T9[. Usually the typical thickness 



of the interface region is of atomic scale, and only the wave function in the outer region 
is of interest. Then, on the quasiclassical level of accuracy, the interface (the knot) can be 
described by the scattering matrix 

In general, the knot may tie together arbitrary number, A^, of ballistic in-channel to the 
same number N of the out-channels. For a specular interface, number of channels N equals 
to 2, and rough interfaces may be modelled by knots with N > 2. 

The waves generated by a source, e.g. on path 1 in Fig.|I], spread to all other paths 1', 2, 
and 2' coupled by the knot. In the presence of an interface, the wave function on trajectory 
remains a valid concept if one interprets the notion of trajectory in a broader sense as a 
set of the points on all the ballistic paths coupled by the knot. For instance, in Fig.|^, one 
understands paths 1, 2, 1', 2' as the parts of a single geometrical object, which we also call a 
"trajectory" . The spatial argument of the wave function will span the generalized trajectory. 
Similar constructions are known in the literature: see, e.g. Ref. where the Schrodinger 
equation is solved on graphs (networks). 

The case of many interfaces requires some preliminary remarks. Consider as an example a 
two layer system Fig.^. //the layers are of the same thickness and the reflections are exactly 
specular, the two outgoing path 1' and 2' meet together again on the upper knot, forming 
a loop i.e. a pair of interfering paths. This causes a major difficulty for the quasiclassical 
theory: Indeed, the envelope function ip obeying the Andreev equation is introduced when 
the phase factor e^^''^^^, C being the distance along the path, is singled out of the full wave 
function. When loops are present and there is more than one path connecting any 2 points, 
the distance £ is ill-defined, and the procedure of constructing the envelope becomes non- 
unique and dubious. Besides, the interference phase factors like e^P^^^i'-^^)/??.^ being the 
lengths of the interfering paths, crucially sensitive to the value of pp and cannot be found 
in the quasiclassical limit where h/pp = 0. 

To overcome the difficulty we note the following: The interference leads to Fabri-Perot 
type geometric resonances and related fluctuations of various physical quantities, perhaps 
locally strong. However, in the limit h 0, the resonances are close to each other in the 
configuration space, and, therefore, the fiuctuations are expected to be effectively averaged 
out when one calculates observables: The latter are given by certain integrals and thus are 
sensitive mainly to coarse-grain features in the configuration space. 

Further, the coarse-grain features (like e.g. the angular-resolved local density of states 
averaged in small volume (S> Xp)) or small interval of directions) are more than likely not 
perceptive to small variations of geometry shifting the positions of the resonances. Hence, it 
seems plausible to assume that the coarse-grain structure can be faithfully reproduced if one 
introduces "virtual roughness", which is small (<C ^o) and not noticeable quasiclassically, 
and performs averaging with respect to the roughness (kind of ergodic hypothesis). In other 
words, on the course-grain level, an ideal surface is expected to be indistinguishable from a 
"virtually rough" i.e. a random surface with roughness W (see Fig.^) small on the typical 
quasiclassical scale, W -C ^o- 

For a rough surface, the picture of trajectories shown in Fig.^ almost never occurs: In 
the quasiclassical approximation, the trajectories are lines with zero (~ Xp) width, and 
the condition that the trajectories 1' and 2' cross each other again exactly at the interface 



6 



(up to ~ Xf), is very restrictive. For this, the surfaces must be strictly parallel and the 
reflections 1'— 3 and 2'— i>4 must be specular (identical) with high precision. Qualitatively, 
the argument here is the same as in the billiard theory where closed orbits are known to be 
rare exceptions. As long as the loops are absent, solutions to the Andreev equation, vary 
smoothly when parameters of the trajectory (e.g. its direction) or the surface roughness 
are changed and have certain limit when the virtual roughness tends to zero. Hence, the 
averaging with respect to the virtual roughness is trivial: it amounts to neglecting it in any 
calculation provided the topology of the trajectories is single-connected. The virtual rough- 
ness (tending to zero) is needed here only as a mean to eliminate the geometric resonances 
which are not of interest because they are not seen on the coarse-grain level of description. 
(Another line of reasoning could be to say that any real sample is always microscopically 
rough so that loops are statistically impossible). 

By these arguments, one comes to the important conclusion that due to the virtual (or 
real) roughness the paths tied together by a knot do not show any further correlations and do 
not (typically) meet each other on other knots. This seems to be an analog to the impurity 
averaging. Effectively, it allows one to average over the Fermi wave length scale from the 
very beginning. 

Uncorrelated multiple collisions with interfaces transform a ballistic trajectory into a 
tree- like geometrical object. To give a general idea of what we mean by a tree, the topo- 
logical structure of one of the possible trees with N=2,3 knots is shown in Fig.^. The tree 
corresponding to a real physical situation will be presented later. 

The main feature of the tree-like trajectory is its one-dimensional character, the property 
which can equivalently formulated as (i) there is no loops or interfering paths; (ii) there is 
only one path connecting any two points of the tree; (iii) the cut of any line produces two 
disconnected pieces. 

Since the tree is effectively 1-dimensional, one is able to repeat Andreev's procedure 
on a tree-like trajectory defining the slowly varying envelope wave ijj{r) by the formula 
^(r) = V'(r)e^P^^(''), where r spans the points on the tree, and C{r) is the coordinate along 
the tree counted off a point. In between knots, the Andreev equation Eq. (|1.6| ) is valid and 
the values of the wave function on a knot are coupled by the scattering S-matrix (see Section 

The purpose of present paper is to extend the existing quasiclassical Green's function 
theory of superconductivity to the case of multi-interface geometry. In essence, the standard 
quasiclassical ("^-integrated") theory of superconductivity is the Green's function version of 
the Andreev equation: Again, the quantum coherence of the electron and hole residing on 
the same trajectory is taken into full consideration whereas the coherence between particles 
occupying different trajectories is neglected. The paths are coupled to each other only by 
the self-consistent effective potentials like various self-energies (impurity, phonon) and the 
pair potential A. The Green function technique has obvious advantages for one is able 
to perform the disorder averaging, include the inelastic scattering and the strong- coupling 
effects etc. 

Although the potential due to crystal imperfections like impurities is far not slowly 
varying, this does not invalidates the quasiclassical scheme if one is interested only in the 
disorder averaged properties. It is well-known that the disorder averaging amounts to the 
impurity self-energy term in the Gor'kov equation which effect is similar to that of the 
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potential energy. The self-energy varies on the same spatial scale as other self-consistent 
potential and as such does not violate classicality. Of course, the imaginary part of the 
self-energy must be small so that the mean free path / is large, I ^ Xp. The quantum 
localization corrections controlled by the parameter h/ppl -C 1 are ignored, which again is 
consistent with the limit h ^ Q oi pp ^ oo accepted in the quasiclassical theory. 

We use the version of the quasiclassical theory |2^j2^ where the main object is the 2-point 
Green's function on classical trajectories. In our opinion, this approach is most adequate 
to the above physical picture of the electron-hole phase coherence spreading along classical 
trajectories. As has already been discussed, in the many-interface geometry the classical 
trajectory becomes tree-like. Accordingly, the arguments of the 2-point Green's function are 
points on a tree. In the present paper, we restrict ourselves to the stationary case, and our 
main concern is the retarded Green's function of the Keldysh technique. 

The paper is organized as follows. In Sect.|l|, we review the quasiclassical theory in the 
formulation based on the 2-point Green's function. The connection to the standard technique 
is discussed in Sect. pi B| . In Sect . [II C| , we briefly show the connection to the Riccati equation 



technique p^ , |25[] , as well as suggest a general method for the case of a periodic potential. In 
Sect.|T|, we derive the boundary conditions for the Green's function on the knot (interface) 
with arbitrary number of channels. In Sect.^, a solution to the multi-channel problem 
of the Andreev reflection as well as the bound states, is given. In Sect.0, we derive the 



interface boundary condition for the Riccati equation. In Sect . [V B| , the boundary condition 
for the Green's function in terms of the transfer matrix is derived. In Sect.0, we show 
the usage of the general approach applying the theory for studying simple examples: (i) a 
film separated by a partially transparent interface from of a bulk material superconductor; 
(ii) two layers of a finite thickness. Motivated by the recent theory of the paramagnetic 



effect I^B], we pay most attention to the case when the phases of the order parameter in 
the two superconductors differ in vr; numerical data for the density of states and superfiuid 
density are presented. The results are summarized in Sect. |Vll| . Details of the calculations 
are collected in the Appendices. In the rest of the paper, h= 1. 



II. TRAJECTORY 2-POINT GREEN'S FUNCTION 

A convenient starting point is the formulation of the quasiclassical technique in terms of 



the 2-point Green's function on classical trajectories; the method was first suggested in |27 



"t-representation"), and in a different form developed in [P^,E3]. The trajectory Green's 



function is introduced via the following representation of the 2x2 matrix Gor'kov Green's 
function lEBI: 



(ri, rs) = -— r-2; ^) + ^ \r,-r2\ ^ ' ^^'""^ " '^sl » 1 

(2.1) 

where mp = pp/v, pp and v being the Fermi momentum and velocity, respectively; e in 
Eq.( p.l| ) is the energy variable (stationary case). For definiteness, we consider the retarded 
Green's function of the Keldysh technique. To simplify notations, we assume a spherical 
Fermi surface; generalization to an anisotropic spectrum is straightforward. 
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Similar to Andreev's procedure, the fast "quantum" oscillations on the scale Xp are 
singled out in Eq. (|2.1|) . Resembling Eg . (11.61 ), the slowly varying quasiclassical envelopes 



g^{ri,r2) obey first order differential equations p8| , p2| , the gradient term of which couples 
only the points on straight lines which are obviously the classical trajectories corresponding 
to a particle on the Fermi surface The trajectory is specified by its direction n and 
arbitrarily chosen initial point R, so that the position r of a point on the trajectory R, n 
can be presented as r = R + xn, x has the meaning of the coordinate on the trajectory. In 
the momentum space, the trajectory n is associated with the points in the vicinity of the 
Fermi surface where the velocity vector is directed towards n. 

For the trajectory specified by {n,R}, one defines the 2-point Green's function 
g^{xi,X2\n,R) H,!!, 



{xi,X2\n,R) 



g+[ri,r2;e) 
g^{ri,r2]e) 



Xi > X2 
Xi < X2 



ri,2 = xi^2n + R 



In many cases we omit R, n and e for brevity and use the notation g^{xi,X2).) 
As shown in p^ , |23|l , the 2-point Green's function obeys the following equations 

d 



IV 



dx-i 



+ H^^{r,) g^{xi,X2\n,R 



g^{xi,X2\n,R) ( -iv 



d_ 

8X2 



ivS{xi — X2) 
iv6{xi — X2) 



r2 



R 
R 



XiTl 



X2TI 



(2.2) 
(2.3) 



where the 2x2 traceless 



matrix H^^, 



H. 



R 



"'£,n ^e,n i 



where A„ is the order parameter (which may dependent on the direction n), and Ps = 

A being the vector potential, and is built of the impurity self-energy and the part of the 

electron-phonon self-energy not included to the self-consistent filed A and . 

The boundary condition to Eqs.( p.2|) , and ( p.3|) is the requirement that g^ is zero at 
\x\ — X2I — ^ C)0, so that g^ is an analytic function of e in the upper half plane for any x\^2 
including \x\ — X2I = 00. 

The advanced Green's function g^ is found from Eqs. (|2.2|) and Eq.( p.3|) with substi- 
tuted for tl^, 

H^ = r,[H^yf, (2.5) 

where is the Pauli matrix and the dagger denotes the Hermitian conjugation. 

Although the observables can be expressed via the quasiclassical 1-point Green's function 
(xi = X2), the 2-point Green's function turns out to be a useful intermediate object. It gives 
a full physical description of the system in the approximation where the part of the orbital 
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degree of freedom is treated classically (no quantum broadening in the plane _L n), with a 
complete quantum treatment of the electron-hole degree of freedom. 

It is important that the construction based on the notion of smooth classical trajectories 
remains valid in the presence of disorder (or phonons), in the standard approximation when 
the scattering is included on the average via the self-energy (provided pi?/ ^ 1, / being the 
mean free path). 



A. Factorization 

To build the Green's function on the trajectory n, R, one first considers solutions to the 
equation 

d - \ 

iv— + H^{x)j(p = (2.6) 

here is a column, ~ (^^^ stands for H^^{r) at the trajectory point r = xn+R. 

Denote if) the row built from a column ip by the following rule: 



Note the identities, 

= -Ma , Ma = , Mb " Ma = (Ma) 1 

By virtue of the identity 

[H^y = -ryH\, (2.7) 
the row 0(x) built from a solution to Eq. (pl6|) , satisfies the conjugated equation 

M) (-^^^ + H^i^)^ = (2.8) 



Combining Eqs.(|2.(j|), and (^.8]), one gets the conservation law, 

^ {Mb) = , (2.9) 

valid for any pair of solutions (f)a{x) and (pb- 

For a general complex e, the Green's function is built of the regular solutions to Eq. (|2.6| ), 
i. e. solutions satisfying the following boundary conditions 

0+(x) ^0,2;^ +00 ,^ 
(f)-{x) ^ , X ^ — oo 

Denote the normalized solutions for which 
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(j)^^\x) (Pf\x) = 1 . (2.11) 

The normalization is possible because the l.h.s. is a (finite) constant as it is seen from 
Eq.(P). 

The Green's function can be written now as 

g [xi,X2)-< (jv) ^-T{N)f ^ ^ [^■^^) 

[ 01 '(Xi) 0^ '{X2) , Xi < X2 . 

Indeed, it satisfies Eg. (p.2|) and Eg. ( |2.3| ) at Xi X2, and is regular at \xi — X2I 00. The 
normalization in Eg. ( P-llD ensures that the discontinuity at Xi = X2, 

g^{x + 0, x) - g^{x - 0, x) = 1 , 



is what is reguired by the 5-function source in Egs. (|2.2| ), and ( p.3|) . 



For clean superconductors with inelastic scattering ignored, S-^*^^) 0, and Eg.( p.6|) 
is nothing but the Andreev eguation Eg. ( p. .61 ). Note that the structure of the eguations is 
not changed when the disorder and inelastic scattering is included via the self energies. In 
this case, however, solutions to Eg. ( p.6| ) have only the meaning of the building block of the 
Green's functions. 



B. 1-point Green's function 

Observables can be expressed via the Green's functions with coinciding spatial arguments, 
and therefore, the 1-point Green's function is the final goal of calculations. 

The 1-point Green's functions defined as g±{x) = g^{x ± 0,x), can expressed via the 
normalized solutions (see Eg. (|2.12|) ) 

g^{x) = 4>f\x)4>i''\x) , g^{x) = 0L'^Hx)0f ^(x) . (2.13) 
This expression can be identically written as 

0_(x)0_|_(x) 0_(x)0+(x) 

where the normalization of the wave functions 0± is arbitrary. 
These matrices are projectors, 

g^g^ = ±g^ , g^g^ = , 9^ ~ 9^ = 1 , Sp^f = ±l (2.15) 

Tagging electron- and hole-like excitations in accordance with the direction of their prop- 
agation (±x directions) and considering examples, e.g. the normal state, one concludes that 
g!^ can be identified as the (guasi)electron part of the Green's function, and g^ is the 
(guasi)hole one (and vice versa for g^). 

Denoting 

a = — , b=—, (2.16 
V- u+ 
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where u± and v± are the components of 0± 
Eg. ( p. 141) becomes 



U±{X) 

v±{x) 
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R 



ah 



. (2.17) 



Another elucidating form of Eq. (|2.14|) is as follows 



where the rotation matrix Oa b 



1 




o-l 



R 



a,b ' 



9- 



O. 



a,b 




-1 



o 



a,b 



1 a 

b 1 



(2.18) 



z.e. 



As discussed in [^,^, the 1-point Green's (".^—integrated") function of the quasiclas- 
sical theory , g^, is given by 

(2.19) 
(2.20) 

(2.21) 



In terms of g^, 



and the relations in Eq.( |2.15D lead to the well-known normalization condition 



m2 



and 



Combining Eqs.( |2.19| ) and (|2.17|) , one gets 
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R 



(2.22) 



1 + ab —2a 
l-ab\2b -{1 + ab) 

This parameterization of the Green's function has been recently suggested by Schopohl 
and Maki |2^ (see, also, [^). The present derivation leads quite naturally to this decom- 
position, and clearly shows the physics behind it. Seeing that a and b may be interpreted as 
the "local" amplitudes of the Andreev reflection for electron and hole (see below) , we call 
them the Andreev amplitudes. 

Finally, the rotation with the matrix Oa,b in Eg. ( p. 181) diagonalizes g^, i.e. 



1 
-1 



6-1 



The advanced Green's function g^ and symmetry relations between g^ and g"^ are dis- 
cussed in Appendix ^ 
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C. Solving the equation of motion 



In this paper we take the approach where the main object of interest is the two component 
"wave functions" (j)±, which factorizes the Green's function and obeys the Andreev-type 
equations. A variety of options can be chosen to find the amphtudes. For future references, 
some of them are discussed in this section. 



1. Riccati equation 

Instead of solving linear equations for two component (j) = ^ ^ ^ ' ^'^^ solves the equation 
for the ratio a = v/u. It follows from Eq. (|2.8|) that a{x) satisfies the Riccati equation, 

ox 

where parameters and A^ are found from the identification 



X 



In the context of the quasiclassical theory, this equation has been first derived by Schopohl 
and Maki Jg^. 

Known a{x), one finds the 2-component function (j){x), 



(j){x) 



const J^^^ ) ( ' / (^"^ + ^''").' ) ^2.24) 



XQ 



To find a±{x) i.e. the solutions to Eq. (p.23|) corresponding to 0± , the Riccati equation 
must be supplemented with the boundary condition which leads to the correct asymptotics 
Eq.fCTD- 



In many cases of interest such e.g. an SNS-structure or isolated Abrikosov's vortex, the 
superconductor is homogeneous at x ^ ±oo. If so, solutions to Eq. (|2.6|) are plane waves in 
the asymptotic region: 



^ conts ^ ^ e^'^""^ 



C,^ = a/ {s^y — A^A*^ , '^^^ > 0. Selecting the waves decaying in the corresponding 
region, one comes to the boundary conditions as follows: 

- 

'■ tt±U=±oo = |x=±00 (2.25) 



An equivalent condition was suggested in [^,0 from "the requirement of the stability of the 
numerical integration procedure" . In the present paper, the boundary condition is deduced. 
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ultimately, from the physical condition that the 2-point Green's function is a regular function 
decaying at large distance from the source. 

The 1-point Green's functions and are found now from Eq.( |2.17l ) and Eq.( p.22[ ) 
with the understanding that 

b{x) = a+{x) , a(x) = 
where a±{x) are the solutions to Eg. ( p. 231) with the boundary conditions in Eq.( p.25|) . 



2. Periodic potential 



In many situations of interest such us vortex lattice, N-S or S-S superlattice, or multiple 
reflections (see below) the potentials are periodic functions of the trajectory coordinate. In 
this case, the Green's functions may be found by the following method. 

A formal solution to Eq. (p.6|), 0(x) = U{x,xo)(p{xo), can be expressed via the evolution 
matrix 



-i J dx' H"(x') 

U{x, xo) = T^e ^0 

where orders the matrices H^{x) in the descending x— order from the left to the right. 

Denote Ul{x) = U{x + L,x) the evolution matrix corresponding to the translation by 
the period of the structure L. As proven in Section 0, the 1-point Green's function can be 
found as 



g'^ix) 



R 



Ul{x) 



Here [ • • • ] stands for the "formating" operation: 



(2.26) 



-R 



Q 



1 



(iR = \l\Q-[ifvQn 



(2.27) 



which returns a normalized traceless matrix |3T[| (similar combination of matrices has 
been introduced in p|). The branch of the square root in must be chosen to satisfy 

9ft {^r{Q\^ > 0. Except for the choice of the branch, Eqs.( p.2^ ), and (|2.271) are same 

for g^. Construction of the evolution matrix Ul{x) in the Riccati equation technique is 
described in Section |B[ 



III. KNOT MATCHING CONDITIONS 



In the quasiclassical picture, particles move on trajectories, usually, straight lines char- 
acterized by the direction of velocity n (and the initial position R). At any point in real 
space, infinite number of trajectories with different n cross each other. Since there is no 
transitions between the intersecting trajectories, the crossings do not lead to any physical 
effect. At some points, called here knots, the quasiclassical condition is violated. At a knot. 
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the particle may leave its original trajectory and continue its motion along a trajectory in 
another direction. In the simplest example of a specular interface Fig. |l|, two trajectories 
1-1' and 2-2' are mixed. In a general case, the knot is a region where transitions between 
in- and out-trajectories are allowed. The in- trajectories (or channels) are those which 
have the direction of the Fermi momentum towards the knot; the momentum direction is 
from the knot in the out-channels (see Fig. ^ |3^. The in- and out-trajectories are somehow 
numbered, / = 1, . . . , A^. We mark by ' the out-going channels so that k' stands for the k-th 
outgoing channels. 

Since the knot is point-like on the quasiclassical scale ~ vp/A, one can talk about the 
knot value of the trajectory "wave function". Denote ipi the 2-component wave function 
on the i-th in-coming trajectory, i = 1, . . . , N at the point where it enters the knot, and 
analogously ipk' i the knot value on the k—th outgoing trajectory. 

The outcome of events happening inside the knot can be generally described by the 
scattering S-matrix. For any specified case, it can be found by solving the Schrodinger 
equation for the electron with the Fermi energy. Here, it is considered as a phenomenological 
input. 

The suggested matching condition reads 



N 

'4'k' = '^Sk'ii)i , (3.1) 

1=1 

where Sku are the elements of the unitary scattering matrix. In the spirit of the quasiclassical 
theory, Sk'i is the normal metal property taken at the Fermi surface; it is an electron-hole 
scalar. This relation generalizes the matching conditions of Ref. to the many channels 
case. 

Taking advantage of unitarity, = S\ the inverse of Eq.(^.lD reads 

N 
k'=l 

Seeing that the conjugated wave function ip always belongs to the second argument of 
the Green's function G(l, 2) ~ ('?/'(l)?/'*(2)), it must obey the matching conditions for ip* i.e. 

N 

^k' = ^k'ii^i (^-^^ 

i=l 

Eq. (|2.6| ) together with the matching conditions in Eq. (|3.1| ), allows one to find the 2- 
component amplitudes on the tree-like trajectory, and, therefore, the Green's functions. We 



remark also that the relation in Eq. (|3.1| ) can be used as the boundary condition to the 
Andreev equation O. 



IV. ANDREEV REFLECTION ON THE KNOT 

In this section, we consider the quantum problem of scattering of ballistic excitations off 
the knot or, in other words, the problem of many-channel combined, Andreev and usual. 
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reflection/transmission. The problem is formulated as follows. On each of the trajectories 
connected by the knot, i, k' = 1,2, . . . N, the order parameter A(x) and, hence, the matrix 
h{x) in Eq.( |2.4| ) is supposed to be known. Since in the ballistic case E = 0, the wave function 
on each of the trajectories satisfies the equation 

iv— + h{x) U = , (4.1) 
ox J 

here x is the coordinate along the corresponding trajectory; this equation differs only in 
notations from the Andreev equation Eq. ( |1.6|) . The scattering of the (quasi)particles off 
the knot is due multiple sequential processes of (i) inter-trajectory transitions described 
by Eq. (|3.1| ), which do not affect the the electron-hole degrees of freedom, followed by (ii) 
intra-trajectory Andreev reflections i.e. rotations in the electron-hole space. The goal is 
to express the amplitudes of the multiple processes via the amplitudes of the elementary 
events. 

On each of the paths, we chose the origin x = at the knot. Then, the coordinate x 
belongs to the region — oo < a: < on the in-coming and to the region < a; < oo on the 
outgoing trajectories. 

First, we consider the plane wave asymptotics at |x| oo where h = const (x). The 
electron- like (hole- like) solution is "i/ei^) = 'ipe^^^^^^ i^hi^) = "^/te"*^^/""), where ipe {'4'h), 
hipe = +^"^6 {hiph = ~^'^h) is the eigenfunction of the matrix h. The eigenvalues ±^ are 
found from i^^l = h^. We supply the energy with an infinitesimal positive imaginary part, 
e e + i6, and impose condition 0=,^ > to specify the branch of \/^. 

The basis for the electron-hole classification is the quasiparticle current 

J,p = i^^f.ij = \u\' - \v\^ (4.2) 

which is a constant of motion, -^jqp{x) = 0, due to the symmetry h'^ = f^/if^. The electron- 
like quasiparticle is identified by jqp > 0. It moves in the direction of increasing x in 
accordance with the sign of the probability current. For the hole-like excitation jgp < 0, and 
it moves towards x = — oo. Note that the solution \|/('^''') are chosen in the way that both 
electron and holes decay in the direction of propagation. 

Below, ■?/'(^'^) denotes the eigenfunctions normalized to the unit flux: 

^(e)t^^^(e) _ I ^ ^Wtf.V^W = -1 , > . 

(The l.h.s. is identically zero in the gap region when < and propagating states are 
absent.) 

Generally, h is x-dependent and the solutions are the plane waves only asymptotically. 
However, the electron- hole classification is unique due to the current conservation in Eq.( |4.2| ). 
One has for the electron- , '^^'^\x), and hole-like, 'i!^'^\x), solutions on out-going (in-coming) 
trajectories 

^^^\x)^{ , . (4.3) 

pie) [ «{e) I ' ^ - U 
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oo or — oo 



^^'^H-)= 1 ^ (4.4) 

where the parameters a*^^''^-' and ^^'^'^^ are found solving Eq.( [4.1|) in the region < x < cx) 
(or —oo < X < 0). 

If considered as a function of x, oiS'^\x) and can be found by solving Eg. ( p.23| ). 

We see that indeed the parameters a[x) and h[x) of the Riccati equation technique have the 
meaning of the instantaneous (local) amplitudes of Andreev reflection, and, therefore, one 
may call them the Andreev amplitudes. 

It generally follows from the current conservation Eq.([4.2|) that 



aW = (a('^))* , |a(-)|2 + |/3(-)|2 = |aW|2 + |^W|2 = l. 
for an open channel, ^ > 0. Seeing that 1/3*-'^^ p = one can enforce 

pie) ^ p{h) 

choosing the overall phase factor in ip^'^''^\ 

The physical meaning of the parameters is clear from Eqs. (|4.3| ), and ( [4.4|) : On the 
outgoing trajectories (0 < x < oo), a*-^-* is the amplitude of the Andreev reflection of the 
(bare) electron injected at x = 0, and /^'•^^ is the corresponding transmission amplitude; 
Q,(^)/^{'*) and are u—,v— components of the quasi-hole having come from x = oo. On 

the in-coming paths, the above is true after the substitution "electron" "hole" . 

Moving towards the knot, quasi-electrons on the in-coming and quasi- holes on the out- 
going trajectories comprise the in-coming states of the scattering problem; the out-going 
states are electrons on the out-going and holes excitations on the in-coming trajectories. 

Let the incoming particle be the quasi-electron approaching the knot along the l-th in- 
trajectory. The source particle generates waves in a// out-going channels. The wave functions 
of the system reads 

vl>(0 = vl>(^) + 5f )vI/W + J2 5f ^ + Yl ^k'^^t' (4.5) 

k^l k' 

where '^^^} and '^'j^^ stands for the trajectory wave functions defined by Eqs.(|4.^), and 



k or k' being the label of the trajectory. The yet unknown amplitudes of the outgoing 
particles, A^, and are to be found from the matching conditions in Eq. ( p.l|) . 

The calculations are most easily done using Eqs.( |C6|) , and (|C7|). It follows by comparing 
Eq. (|C6|) with Eqs.([4.3|), ([4.4|), that one may put Um^^i = ain and fik' = ai*^^- From Eq. ( |C6|) 
and Eq . (|C^ ) , one sees that the wave functions on the source trajectory at x; = must be 
proportional to ( le)) , where 



= (/|Std(^)5,|0. (4.6) 



Here and below, Si, is the full S-matrix taking into account multiple events of the Andreev 
reflection. From Eq.( 
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Si = (^S^ - ' ; (4.7) 

where the diagonal matrix with the elements (d*^^'''^) kk = a^i^''^'' and superscripts 

means the //-element must be put to zero. By (/|Z|m), /,m = 1, ... we denote the 
matrix element Zim- 

The parameter Oq^^ in Eq.( [4.6|) has the meaning of the amplitude of the Andreev backscat- 
tering of a bare electron by the knot as a whole. From the condition that the wave function 
has the u — v structure at x; = like ( fe)), one finds Bi, i.e. the amplitude of the Andreev 
reflection of the incident electron excitation. After some algebra 



Here, the denominator can be understood as due to multiple Andreev reflections p3 |. 

The wave function + Bf^'$['^^ at = equals now to C(^Je)) where C = 

pi^^* (l - . Looking at Eq.(^, one find the rest of the scattering amplitudes: 

Ai} = - \--^{k'0(^^SiP^^>\l) , = - j--^(A;|/3W^t^(^)5,/3('^)*|/) . (4.9) 

1 - a'Qi'a] 1 - a'o/a] 

Similarly, one derives the scattering amplitudes for the quasi-hole coming to the knot on 
the n'-trajectory. Analogously to Eq.( |4.5|) , the wave function. 



k'^n' k 

contains the scattering amplitudes which are found from the matching conditions. The 
corresponding expressions can be obtained by the substitutions: (e) ^ (h), I ^ n', and 
5" S^, and Si S^,, 



sl = (s- (d(^))("')^d('^) 



For the hole incident on the n'-trajectory, the amplitudes of the Andreev reflection, \ 
scattering to the hole state on the k-th trajectory, A^^ \ and scattering to the electron state 
on the k'-th trajectory, 5^? read respectively 



B' ' = 

n' 


1 


(h) (e) 






1 


1 


(h) (e) 






1 


1 


(h) (e) 



(4.10a) 

{k\(3^''^sl(3^''>\n') (4.10b) 
{k'0^'^Sa^''^sl^^^>\n') (4.10c) 
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The presented formulae give the amphtude of scattering from a propagating channel 
to another propagating channel. The scattering of the excitations is a result of multiple 
sequential events of two types: (i) on the knot inter-trajectory transitions described by the S- 
matrix in Eq. (|3.1|) , and (ii) intra-trajectory processes of the Andreev reflection/transmission 
with the amplitudes a'^^'^'' / (3^^'^\ Expanding the effective S'-matrix Si in Eq.( [4.7[ ), Si = 
5 + 5(aW)«5td(^)^ + 5(dW)(')5td(^)5(dW)«^td(^)5 + . . . one sees that the full amplitude 
of the scattering event n' <— m is the superposition of all different paths connecting the 
initial and final states with electron hole transformation on each step. 

The theory gives exact amplitudes of the multiple scattering expressed via the amplitudes 
of the elementary processes: the normal metal S-matrix and the intra-trajectory Andreev 
amplitudes. In the simplest case, when N = 2, and A = on two out of the 4 trajectories, 
the above formula reproduce results of the theory of Andreev reflection in the NIS structure 
21. 



A. Bound states 

Bound states are physical solutions existing in the absence of a source. The physical 
solutions are those when the matching conditions on the knot are simultaneously satisfied 
with the requirement that the wave functions decay far away from the knot. The electron 
and holes states defined earlier (with > 0) have the property that they decay in the 
direction of their propagation. Therefore, the wave function of a bound state \l/bound has the 
form 

Abound = J^i^fe^f^ + J2^k''^i' 

k k' 

where the coefficients A's and 5's are found from the matching conditions. Again, looking 
at Eqs. (|4.3| ), and ( [4.4|) one sees that in Eq. (|Cl|) , Hk' niay be identified with Q!^^\ and z/j with 
Then, Eq. 



{a^'^}) =det 1 - Sa^'^^S^a^'^ = , (4.11) 

gives the condition for the wave functions to be matched on the knot. The Andreev am- 
plitudes a^^'^ and a^'^^ are functions of energy e, and the bound states exist at the energies 
where Eq.( [4.11|) is satisfied. 



B. Example: Rough surface, anisotropic superconductor 

The rough surface reflects waves in many direction. As the simplest model, we assume 
that the surface reflection couples together only 2 in-coming directions " 1" and "2" to two 
outgoing "1"' and "2"'. The model corresponds to a = 2 knot. In what follows we 
calculate the amphtude of Andreev reflection by the knot and consider the bound levels. 

The unitary 2x2 scattering matrix of the knot may be taken in the form 
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provided -Ri + -R2 = 1, where Ri = |rip {R2 = |?"2p) is the probabihty of reflection 1-^1' 
(2-1')- 

Given the profile of the order parameter, one can find the wave functions, and the 
Andreev amphtudes a^^''^^ and Here, the matrices 

are taken as input, each of the a's is a functions of energy. 

The energies of bound states are found from Eg. ( |4.1l| ), which takes the following form 



V{e) = - aS^^aS?^)(l - ^aj^) + i?2(l - aPa'^^^){l - a^'^^aS?^) = . (4.12) 

The bound states exist only in the gap region at the energy interval where \a^^l \ = 2' I ~ 

Essential physics can be grasped by the simplest model where the order parameter A„ 
is a constant at each of the trajectories: A„ = Ae"^". Then, 



where -0^ is a function of energy, e'''^^ = — Ap — e^) /{e + A^ — e^). 
Eq. (|4.12| ) is conveniently transformed to the form, 

COS[lpe + ipw + ip22') = Rl COS + R2 COS , (4.13) 



One sees that the existence and position of the bound state is sensitive to the surface 
roughness only if either the incoming or outgoing channels are not equivalent i.e. ipu = 
'Pi — V2 ^ 0, or {pi/2> = v^i' — ip2' 7^ 0. In other words, mixing of identical channel does not 
affects the levels. 

Consider now the possibility, which may exist in the case of a d-wave superconductor, 
that the order parameter changes its sign on the 1^1' and 2^2' trajectories. A smooth 
surface mixes only trajectories with close transverse momenta; then the trajectories are 
almost equivalent and their coupling does not shift the levels. On the contrary, a backward- 
like scattering splits the degenerate levels: In the model under consideration, the backward- 
like scattering corresponds to the phase factors (pi = (^2' = and !f2 = 'Pi' = 0. Then, from 
Eq.( [4.13| ) cosipe = R2 — Ri- The bound state energies are 



Abound = ±/R2 A . (4.14) 

One concludes that the presence of substantial spectral weight at low energies is not likely 
if scattering in the backward directions is present: ~10% probability the scattering moves 
the levels from zero energy to ~ 0.3A, of the order of the gap. 

The amplitudes of scattering of excitations can be found from Eqs.(|4.8|), ( |4.9| ) and ( [4.10c| ) 



As an example, the amplitude of the Andreev reflection of the electron-like excitation inci- 
dent on the trajectory "1", b[^\ reads 
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where the following notations are used 



~ih) «1 — <^2 

"1 1 (e) {h) 

1 — ayai 

(e) (e) 

(e) ^ Qj/ -^2^ 

"l' 1 (e) {h) 

(e) (e) 



a 



(e) (h) 



1 



The shortest way to derive this result is to apply the rotation transforming and a^i to 
zero as explained in Section 0. 



V. MATCHING GREEN'S FUNCTIONS 



As has been discussed in Section [II A| and |1I B|, the Green's functions can be built from 



the regular solutions to the Andreev equation Eq. (|2.6| ). When the trajectory coordinate 
X extends from —oo to oo, the regularity requirement leads to the boundary conditions in 
Eq.( p.lO ). In the case of a trajectory ending in or originating from a knot the boundary 



conditions must be reformulated. 

First consider an isolated knot mixing semi-infinite trajectories (with no more knots on 
them). With the origin chosen at the knot, the trajectory coordinate x„ extends from — oo 
to on the ra-th incoming trajectory, and < x^' < oo on the /c'-outgoing one. As before, 
the requirement, 

(/)_,m(-oo) = , </)+,fe/(oo) = ,m. A; = 1, . . . , , (5.1) 

uniquely (up to a normalization factor) defines the solutions 0_^„(x„) and (pj^^yi^k')- Denote 
the knot values of the regular solutions as 

0-,m(a;m = 0) = ( ^™ ) ' = 0) = ( , ) ,m, A; = 1,...,A^. (5.2) 



. 1 / V^k>, 

For convenience, the normalization is chosen so that one of the components equals to 1 at 
the knot; the parameters or h^i are "bulk" properties independent on the knot. 
The problem in hand is to find the knot values 

0+,,(xi = 0) = , 0_,„/(x„/ = 0) = ' l,n=\,...,N 



which give the boundary condition to Eq. (|2.6|) needed to evaluate (f)+,i{xi < 0) and 

0_,„'(X„/ > 0). 

To find 0+^/(0), one notes that by virtue of the matching conditions in Eq. ( |3.1| ) and 
Eq.( |3.2|) , a finite 0+^;(O) generates waves in all other channels, outgoing and incoming. In 
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a regular solution, all the secondary waves must decay while propagating from the knot. 
This condition fixes the u — v structure of the secondary waves: in each of the channel, the 
incoming m 7^ / and any outgoing one fc', the generated 2-component wave functions (at 
x = 0) must be proportional to that in Eq. (|5.2| ). As proven in Section 0, the matching 
condition allows one to find the u — v structure in one of the channels provided, as is the 
case here, it is known for all other channels. 

Changing notions in formulae in Section y (/ifc' hk' , i^i^i —>■ ai , ui-i = bi), one gets 
from Eq.([U4D 



bi = {l\S^bSi\l) (5.3) 

where 



1 



a = diag(ai, a2, . . .) and b = diag(6i/, 62', • • •)) the superscript has the meaning that the 
l-th element on the diagonal must be put to zero; and (/|(. . .)|/) = (. . 

Repeating the arguments, one finds the boundary value a„/. Changing notations in 
Eq.( |C8|) (/i^/ ^ an'), one gets 

o-n' = {n'\SaS^^,\n') . (5.4) 



where 



From the derivation in Section 0, it is clear that both Eq.( |5.3| ) and |5.4| are just different 
forms of Eq. (|C3|) , which reads in the present notations 

V{{a}, {b}) =det 1 - SaS^b = . (5.5) 

This equation should be understood in the following sense: 

Suppose one seeks for the boundary value of 6; for the /-the in-channel. Then, one 
formally solves Eq. ( |5.5|) relative to ai, the obtained value (the inverse of the r.h.s. Eq. (|5.3|) ) 
gives In the same manner, one finds the knot value of a/.' on the /c'-outgoing trajectory as 
the inverse of the root of Eq.( |5.5| ) relative to bk'- The procedure does not pose calculational 
problems since the determinant is a linear function of any of a's or 6's. Eq. (|5.5| ) represents 
most concise and symmetric form of the boundary condition to Eq. ( |2.6| ) . 

Summarizing, the Green's functions on trajectories linked by a knot is calculated in the 
following scheme. First, one solves Eq.( p.6| ) with boundary condition in Eq.( |5.1| ) on each of 



the trajectories and calculates functions (f)^^rn{x < 0) and (j)+^k'{x > 0); the parameters 



and bk' in Eq. ( [5.2[ ) are then also known. The next step is to calculate the knot value of 6's on 
the incoming trajectories and a's on the outgoing ones. This is done by formulae in Eq.( ^.3|) 
and Eq.(|5.4|). Having obtained the boundary values, one solves Eq. (^.6]) for (x < 0) 
on the incoming trajectories and 0_^fc'(x > 0) on the outgoing ones. The 1-point Green's 
function is then built from (j)± by the recipe in Eq.( |2.20D . 



22 



In the Riccati equation technique, one first finds the Andreev amphtudes am{x) and 
bk'{x), m,k = 1, . . . ,N from Eqs. (|2.23| , |2.25|) . Then, Eqs. (|5.3| ), and (|5.4| ) provide the initial 



value for hmix) and a^^ix), solutions to the Riccati equation. The Green's function is then 
given by Eq.(p:22D. 



The matching conditions can be also expressed via the transfer matrix as derived in 
Section [VB| and in the case of a = 2 knot explained in detail in Section ^ 



This scheme is also applicable when the trajectories connected by the knot under con- 
sideration may enter other knots. As a matter of principle, one assumes that the system 
under consideration is finite, and it is surrounded by a "clean" material where trajectories 
are infinite lines without knots. Then, one solves the problem for the knots on the boundary 
and moves inwards towards the knot of interest. In the one-dimensional topology of the tree 
with only one path connecting any two knots, the procedure is unique. 

A. 2 X 2 case 

The most simple case is when the knot mixes two incoming and to two outgoing trajec- 
tories (A^ = 2) as e.g. in case of specular refiection on an interface. The unitary S-matrix 
coupling 1 and 2 incoming trajectories to V and 2' outgoing ones (see Fig.^, may be taken 
in the form 

r s \ 1 i2 _ D I „|2 



S=\^_^*^*j Ar\' = R, \s\' = T . R + T = l. (5.6) 

Here, r and s are the amplitude of the process 1 ^ 1' 2 ^ 1', respectively. 
Presenting the wave function on each of the trajectories at the knot as 

«1(2) ^ ( 1 



the matching condition in Eq. (|5.5|) gives the following relation between the parameters 

R{1 - aibr){l - a2b2') + T(l - ai62')(l " a2&i') = , (5.7) 

which serves as the boundary condition for Eq. ( p.6| ) or Riccati equation Eq.( p.23D . 

The usage of it has been explained in Section ^ Reiterating, the parameters ai,2 {bi',2') in 
Eq.( |5.7| ) are found from the regular solutions to Eq.( |2.6| ) or Eq.( |2.23| ). They are independent 
from each other and the properties of the knot. The actual meaning of Eq. (|5.7| ) is that when 
it is resolved relative to ai^2 (&i',2') the inverse value gives the initial condition 61,2(3; = 0) 
(ai',2'(0)) I.e. 

, R{1 - a2b2')bi' + T(l - a2bi,)b2' , . 

''^'^ = i?(l-a2M+T(l-a2M ' 

R{1 - a2b2')ai +T{1 - aib2')a2 ^, 
<"> = R(l-a,b,)+T(l-aM ' <^-'" 

and the expressions for 62 and 02' obtained by the substitution 1^2. 
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B. Transfer matrix 



Sometimes it is convenient to consider a pair of trajectories, tag them to 1 and 1', as 
pieces of a single trajectory (see Fig.^). We assign a; < to the path 1 and a; > to 1'. 
Then Eg. (|2.6| ) is vahd for any x excepting the knot point x = 0. The knot at the trajectory 
1' ^ 1 is included via the 2x2 transfer matrix Aii>^i: 

0(x = +0) = Mi'^i(pix = -0) , ^{x = +0) = ^{x = -0)M^,l_^ (5.10) 

as explained in detail in Appendix Sect. |^. The transfer matrix is found from the require- 
ments that (i) the matching conditions in Eg. (|3.1| ) are satisfied; (ii) waves on the trajectories 
other than 1 and 1' are regular. 

Denote (p+{x > 0) (0-(a; < 0)) the solution to Eg. ( p.6|) regular at +oo (— oo) as in 
Eg.( |2.10D . The transfer matrix allows one to continue the solutions across the knot: 

0_(+O) = 0_(-O)M-,Li , M~0)=MyUM+0) (5.11) 

In accordance with Egs.( |2.2lD , and ( |2.14| ), the 1-point Green's function on the trajectory 
1' at the knot can be found as 

1 ^ ,R^ _ 0+(+O)0-(+O) 
2^ +'"^"04+O)0+(+O)- 

Applying Eg.( p.llj ), one gets from here that 

hi , s«i = ^+(+o)0-(-o)-Mr/^i 

2^ ^''^ 0_(_O)A1-,L,0+(+O) ■ 
Similarly, for the trajectory 1 

1 ^ ■Mr,ii0+(+o)0-(-o) 

2^ ^_(-0)M^,l^M+0) 
From here g^Air^i = A^i'^i^f or 

For an arbitrary interface, this relation gives the boundary condition for the quasiclassical 
1-point Green's function. With the help of Eg. ([El|) or Eg.(^), the transfer matrix A4 is 
generally expressed via the Green's function on the other trajectories coupled by the knot. 
In the next Section |V B 1| , we the explicit expression for the transfer matrix is presented for 
the simplest case of 2 in- and 2-out channels. 



1. 2x2 case 

In the most important case of a 2 x 2 knot Fig.|^ (e.g. a specular interface), the transfer 
matrix can be found usual the general formula derived in Sect.^ A more simple way is to 
make the derivation from the scratch in a specially selected basis (see Sect. [EOal , for details). 
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For the knot with the S-matrix in Eq. (|5.6|) , the transfer matrix Eq.( [E5|) and its inverse 
read 

a^.-.-^(i-t|i^. (5.13) 



where g2',2 is the normahzed (^|/,2 = 1) "across-knot" Green's function. It can be presented 
is different forms. 

Its matrix structure is most transparent when ^2'»2 is written in a factorized form as 

^ (1 + ^2'.2) = ^</'2',+ 02- , N = 02,-<^2',+ (5.15) 

where 02',+ and 02,- are the functions introduced in Sect . |II A] taken at the point adjacent to 
the knot on the trajectory 2' or 2. They do not dependent on the knot parameters R and 
T. One may think of ^2'»2 as a 1-point Green's function on the virtual trajectory built of 
the pieces 2 and 2'. 

From Eqs.( ^.15D and (|2.13| ), one concludes that (1 + ^2'«2) oc (1 + g^){l + ^^) 



Eq. ( |5.15| ) can written in terms of the Andreev amplitudes Eq. ( |2.16| ) , as 



1 / 1 + 0-2^2' -2a2 



or 

= [ 262' -(1 + a262') J • ^^-^^^ 

The "across interface" Green's function can also be written as 

92',2 = ^ , — (^2' +^2 + J[^2',^2]- ) , (5.17) 

1 + 2i^2',^2j+ V 2 y 

where ^2,2' are the knot values of the usual 1-point Green's function on the trajectory 2 and 
2'. One should realize that unlike 0± in Eq. (|5.15| ) and a, V in Eq.( ^.16| ), both ^2 and ^2' are 



modified by the knot scattering, and only their combination ^2'»2 is knot independent. 

Using the transfer matrix approach, one can derive the boundary condition to the Riccati 
equation on the N=2 knot. Most easily this can be done using the transfer matrix in Eq. (|E4|) . 
Same the result one can get from Eq. ( |5.7| ). 

We have just presented the boundary condition for the Green's function on an inter- 
face which mixes 2 in-coming and 2 out-going trajectories (e.g. for a specular interface): 
the Green's functions on the interface are linearly related by Eq. (|5.12| ) (and the analogous 
relation for the channel 2 and 2') where the transfer matrix M. and M.^^ can be found 
from Eqs. (|5.13|) , (|5.14|) and ( ^.1?1 ). Using these relations, one is able to re-derive Zaitsev's 



boundary conditions H for a specular reflecting interface. 
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VI. MULTILAYER SYSTEMS 



The purpose of this section is to show the usage of the general theory in practical calcu- 
lations. First we consider simplest geometry that is a layer deposited on the flat surface of 
a bulk material with a partially transparent interface. Together with the totally reflecting 
outer surface, there are two coherently reflecting planes. The other geometry is a system of 
two layers of arbitrary thickness in contact, in which there are three reflecting planes and 
rather complicated picture of multiple scattering. 

Since out main intention is demonstrate how to use the general formula, we allow our- 
selves not to worry about the self-consistency of the pair potential. For simplicity, we 
consider the ballistic case Sjmp = 0, and the pair potentials in the left (1) and right (r) 
regions are taken constants A/ and A^. 

A. A film 

The tree-like trajectory near the interface between a layer of thickness dr and semi-infinite 
space is shown in Fig.|^(a). To build the tree, one considers a particle coming along the path 
(at the angle 6) marked in Fig.0 by "1" which denote both the location and direction. Due 
to the partial reflection, a wave on the trajectory "4" is generated. The waves on the paths 
"2" and "3" are generated due to transmission. The paths "2" and "3" are the semi-infinite, 
whereas the trajectory "4" comes to the interface again as "5" (the total reflection does not 
interrupt motion in between "4" and "5"). Again, waves on "6" and "7" are generated, and 
the path continues towards "9" etc.. The topological structure of the tree-like trajectory is 
presented in Fig.|^(b). 

To find two-point Green's function g^{xi,X2), one solves Eq.( |2.3| ) where the coordinates 
Xi^2 correspond now to the points on the tree Fig. 0(b) with the understanding that the tree 
coordinate x includes information about both the position and direction of the momentum. 
Due to the one-dimensional topology of the tree, the method described in Sect.|I| is directly 
applicable. As before, the 1-point Green's function g^{x) is given by Eq. (|2.19| ). 

The matrix in Eq.( p.3| ) is either 

for the tree coordinate x in the left or right regions. For future references, the free bulk 
1-point Green's function in the left (right) region equals 

— ^ iiR 

where = ,/{^TW^W^, ^zf.) > 0. 

Considered as a function of Xi, g^{xi,X2) has a source at Xi = X2 which generates waves 
propagating away from X2. The regularity condition requires that the waves decay when 
propagating from the source to branches of the tree. The propagation in between the knots 
is described by Eqs.( p.6|) , or ( p.8|) , and the knots are incorporated by the matching conditions 
in Eq.( p.lD or their more advanced version in Eqs. (|5.5| ), or (|5.10|) . 
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Let us first find 1-point Green's function at tlie tree point x in between "4" and "5" . In 
accordance witli Sect.|l|, one lias to find solutions 0+ which describes the wave spreading 
from the point x in the positive direction, and 0_ propagating in the opposite direction. In 
the present example, the wave 0+ spreads to the paths "5", "6", "7", "8". . . , and 0_ spreads 
to "4","3","2","1". . . . We chose to think that the particle moves along the "root" path 
"1"_^"4"_^"5"_, "8"^"9"..., and exclude the "side" branches "2","3", "6" ,"7", ...using 
the transfer matrix approach (see Eq.(|5.11|)). 



Take e.g. the knot where the trajectories "1-4" meet (see Fig.^(a)). The transfer matrix 
Ma ^1 can be expressed in accordance with Eq.( |5.13| ) via the "across-knot" Green's function 
g^,2- In the present simple case, when "2" and "'3" extend to infinity and is same for 
"2" and "3" , one can conclude from Eq. ( ^.15|) or Eq. (|5.16|) that g^,2 = Qo i , where g^i is the 



bulk Green's function in the left region. For any of the identical knots, the transfer matrix 
reads 



M = ^ — 1 g^j 

where R = |rp and T = 1 — R are the interface refiection and transmission probabilities. 

The functions (f)±{x) on the root trajectory, where x is the coordinated along the root 
counted from a knot, are found with the help of Eq. ( p.6|) supplemented with the boundary 
condition connecting the 2-component amplitude leaving the knot (pout (out = "4" 
"12", . . . ) via the incoming wave (in = "1", "5" ,"9", . . . ) 



■'out 



Mcpi 



In the present case, when the free motion on the root trajectory is perturbed by the equidis- 
tant knots, one can use the method developed in Sect.^ for periodic potentials. The period 
of the structure is 2De, Dq = dr/ cos 9 where 9 is the angle between the direction of the 
momentum and the perpendicular to the interface. 

The functions 0±(x) are eigenfunctions of the evolution operator U2Dg{x) generating 
the translation by the period x ^ x + 2D0 (see Sect.^). The free evolution operator 
U^'^\x + xo,xo) in the right region is 

U^''\x + Xo,Xo) = e-^^^"" = cos(^) + i^J,,sin(^), 

^1^^ being the bulk Green's function. 

The full evolution operator U2Dg{x) reads 



Uo{x) = Aexp (^^0^;.^) (1 - {^^ali^De 



where A = Finding the two eigenfunctions of this matrix, one knows 4'±{x) and, 

therefore, the full 2-point Green's function from Eq.( |2.12D . 



As explained in Sect.|B[ the 1-point Green's function can be extracted from U2Dg{x) by 
purely algebraic transformations. The Green's function for the direction of the momentum 
{p)z = Pf COS 9 at the distance from the interface z {z > in the right region) reads from 
Eq.(Pl 



27 



R 



R 



exp X 



V 




X 



(6.1) 



where the "formatting" operation [ . . . ] is defined in Eg. ( ^.271 ) . The "formatting" can be 
performed analytically but the result looks rather awkward and hardly any information can 
be extracted from it without a computer. On the other hand, the "formatting" operation 
is easily implemented numerically, and for this reason we leave as final the expression for 
Green's function in Eg. ( |S.1D . 

Consider now the left region and the knot "1"- "4" in Fig.|^(a). The left region Green's 
functions are those on trajectories "2" and "3". To apply formula in Sect. |^ and one should 
substitute 1 for "2" and 1' for "3". Since trajectories "2" and "3" are semi-infinite, the 
combination 0+(+O)0_(— 0) is proportional to the bulk value (1 + gQi)- The transfer matrix 
A^4^i contains the across- knot Green's function ^4,1 analogously to Eg.( |5.17| ). It is easy to 
see that ^4,1 eguals to just found g^{z = +0, 9). Therefore, the Green's function on the left 
side of the interface is 



-0,1 



-R 



1 



T 



1 + R 



g\z = +0, 



At other points in the left region (2; < 0), the Green's function is found with the help of the 
free evolution operator. 



exp ( '^-^gQjX ) g^[z 



-0, 6) exp 



(6.2) 



x=\z/ COS i 



i.e. $5 g^{z 



In Fig.|], we show the density of states on the film side of the interface, 
+0, 9) Eg.( |6.1| ), for Dq = v/\Ai\ and the pair potential in the left and right parts of different 
signs, = — A.^; the curves parameters differ in the reflectivity R increasing from zero in 
Fig.|(a) to R=0.9 in Fig.|(d). 

When R = 0, one sees in Fig.§(a) two (zero width) peaks in the gap region \e\ < |A|. 
The peaks are due to the bound states well-known known in the theory of anisotropic 
superconductors ||3^ (see also Sect . [IV B]) . The e = bound states exist near the trajectory 
point where the phase of A changes abruptly by vr. When the thickness dr is flnite, the 
levels are at a flnite energy p6[ due to the overlap of the wave functions (e.g. of the states 
on the "2"- "4" and "5"- "7" paths in Fig. 0(a)) and the level repulsion. The overlap of the 
separated in space levels and, therefore, the level splitting are exponentially small when Dg 
is large. 

When R is flnite, the splitting increases. First, the reflection gives rise to the on-knot 
overlap of the levels belonging to the the same knot, e.g. the "2"- "4" and "1"- "3" levels". 
By this mechanism, the level is split to ±a/R|A| (cf. Eg.( 4.14] )). Second, the on-knot 
overlap in combination with the next neighbour overlap discussed earlier, mixes together all 
the bound states and transforms the discreet levels into bands. This behaviour is clearly 
seen in Fig.||(b)-(d). 
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B. Sandwich 



In this Section we consider a more general case when the left region is a finite layer of 
thickness di. As previously, the order parameter is assumed to be constant in the layers. 

The typical tree-like trajectory formed by multiple refiections on the outer surfaces and 
the interface, is shown in Fig.^(a). As in Fig.^ the numbers tag the coordinate on the 
trajectory. Topological structure of (a fragment of) the tree is shown in Fig.p(b); the tagging 
in same is in Fig.^(a). The centre of the tree is (arbitrarily) chosen at the "5"- "8" knot; 
the tree structure looks same if viewed from different knots. The pieces of the tree with the 
arrows in the horizontal direction correspond to the the left layer, and points on the vertical 
lines belong to the right layer. Generally, the tree-like trajectory covers (almost) all space 
but remains nevertheless topologically one-dimensional: The features discussed before are 
clearly seen here that is (i) if a line of the tree is cut, two disconnected pieces are produced 
or, equivalently, (ii) there is no closed loops on the tree. 

First we calculate the knot values of the Green's functions, for the central knot "5"- "8". 
Other knots are equivalent to the central knot. On both horizontal and vertical branches 
in Fig. 1^, the arrays of knots are periodical, separated by 2Difi, Di e = di/ cos 9 for the 
horizontal branches (the left layer) and 2Dr^e, D^fi = dr/ cos6, on the vertical branches (the 
right layer). 

As in the previous section (see Sect. |B| for prove), the 1-point Green's function at "5", 
g^, is simply related to the evolution operator advancing the wave function at "5" to 
the periodically equivalent point "9" (see Fig. P)(b). Crossing the knot from "5" to "8" 
with the help of the transfer matrix, build analogously to Eq. (|5.13|) , 



2r* V 1 + 

and moving from "8" to "9" by exp{2iDg rgor)^ o^i^ get as the ordered product of the 
two matrices. The same matrices but multiplied in the different order, give the evolution 
operator f/s^s and, therefore g^. 

Changing notation in Eq.( p.l7| ) and collecting formulae together, one gets 

97.6= ^ , ir^ ^1 (^7 + ^6 + ^[^7,^6]- j (6.3a) 
1 + 2197,961+ V 2 y 



1^ = J=. 



R 



R 



1 - YT^^^'O ^^P*'^^^^''^^^/^'* 



(6.3b) 
(6.3c) 



These equations allow one to find the knot values of the Green's function in the right region 
via the left region counterparts. 

In the same way one can derive expressions where g^j are related to g^^^ 

98,5 = —TT^-r^ \ 9s + 95 + hgs, 95]-] (6.4a) 
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9^ 



9e 



-R 



T 



.6 exp{2iDg^igQi/i 



R 



(6.4b) 
(6.4c) 



Eqs.( |6.3| ) and ( |6.4| ) allow one to find iteratively the knot values of the Green's function. 
Unless the reflection R is too small, the iterations converge rather fast. For almost trans- 
parent interfaces, i? -C 1, a slightly different procedure is more efficient: as the periods, one 
chooses the paths like "4" ^"5" ^ "7" ^ "14" ^ "16". 

Given the knot values, the Green's function at other points can be calculated by formulae 
analogous to Eq.(|6!^ 



Fig.|^ shows the trajectory resolved density of states at the interface, 9 g for the 
case when the A; = — A.^ and the layers of equal thickness Di e = Dj.q = v/\Ai\. 

As expected, the sandwich with a transparent interface, R = 0, has a considerable 
spectral weight at low energies which is represented by the band centred at e = (see 
Fig.|lO|(a)). The overall picture is very different from the BCS density of states: the spectrum 
is given by well-defined bands with strong edge singularities. As in case of a film, the 
reflection splits the e = bound states, and the bands move towards higher energies. When 
the reflectivity is as low as 0.1 (see Fig.|TD|(b)), there is no states at, and in the vicinity 
of e = 0. The forbidden bands become more narrow, and the edge singularities become 
smoother. From Fig.|TD|(c) and (d), one sees that for R > 0.5 the states are pushed to the 
energies > A. 

In the next section, we use these results to evaluate the "superfiuid density" , an observ- 
able sensitive to the shape of the density of states. 



1. Superfiuid density 

In this section we calculate ps, a parameter which controls the current density j induced 
by a weak spatially homogeneous static vector potential. A, 



c 1 



Xl being the (bulk) London penetration depth at zero temperature. In the two-fluid lexicon, 
Ps is the "superfiuid density" or the "fraction of superconducting electrons" . 

In the present case of a two-layer system, the local current induced by in-plane homoge- 
neous vector potential is 2;-dependent, being proportional to the local density of states. The 
total current through the layers is proportional to the average, 

dr OO 

dz ps{z) , ps{z) = l- I de ( u{e,z) (6.5) 



di + dr J J \ de 

-di -OO 

where /o is the Fermi function, and z/(£, z) is the local density of states, 

iy{e,z) = Re J {g^{e,n, z))^^ . 
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The averaged superfluid density ps in Eq.( |6.5|) is conveniently written as 

Of) 

where i^i^ri^) is the averaged density of states in the left (/) and right (r) layers. 

To calculate z/; r(^), one finds the Green's function, as explained in the previous section, 
and perform integrations with respect to the coordinate z and the direction n. The spatial 
dependence, found from the knot values by formulae analogous to Eq.( |6.2| ), is simple and the 
z— integration can be done analytically. The averaged density of states in the right region 
reads 



i 



+ '-^9oA9o,r, 9^]- + ^ (^-^^) i9o,r, 9^]^)^ (6.7) 



where, = cos 9, 7^,^ = 2dr^^/ /if, Qq^^ is the bulk Green's function, and is the knot value 
of the Green's function on the tree corresponding to the angle 6 (see previous Section). After 
the substitution, r — > / and — > g^, Eq.( |6.7| ) gives ui. 

The integration with respect to fi in Eq. (|6.71 ) and e in Eq.( |6.6D can be performed only 
numerically. The integration in Eq. ( |6.6| ) along the real e axes may be slowly converging 
due to the band edge singularities; for better convergence, one may integrate along line 
Q e = or transform the integral to the Matsubara sum. 

We evaluated numerically the superfluid density for a sandwich with equal thickness of 
the layers di^r = v/\Ai^r\ and the vr differences in the order parameter phase = — A^. In 
Fig. 11, the superfluid density as a function of temperature is shown for different reflectivity 
R. 

The curve for R = shows large negative ps at low temperatures which would lead to 
amplification of the applied magnetic field rather than the Meissner screening. This feature 
is due to the large low energy spectral weight seen in Fig.|l^(a). Therefore, our data support 
the recent idea put forward by Fauchere, Belzig, and Blatter |^ about the paramagnetic 
instability near the surface where the order parameter changes its sign. However, one sees 
m Fig.0 that the effect is very sensitive to the presence of the partially reflective interface: 
reflection with the probability as low as 4 percent makes ps positive at any temperature. 



VII. CONCLUSIONS 



In this paper we have reconsidered the part of the quasiclassical theory of supercon- 
ductivity which concerns interfaces between superconductors (SIS) or a normal metal and 
a superconductors (NIS). Since the interface violates the condition of applicability of the 
quasiclassical approximation, the reflection and transmission processes must be included via 
a boundary condition. In the approach taken in the paper, the master boundary condition 
in Eq. (|3.1|) is formulated for the effective wave functions factorizing the 2-point Green's 
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function. In the boundary condition, the two-component amphtudes in N in-coming and 
N out-going channels are related to each other via the S-matrix. The latter is sensitive 
to microscopic details of the interface and is considered as an input in the quasiclassical 
theory. The theory is equally applicable to specular interfaces {N = 2), as well as to the 
many channel case which models a rough surface or interface. In Sections |IV|, |V|, and |V B 



the master boundary condition is reformulated in various forms, suitable for the one or the 
other application. 

In Sect.^, we have presented a general solution to the ballistic problem of the scattering 
of electron-and hole-like excitation. This result extends the theory of the NIS interface |^D[ to 
the many channel situation; SIS case is also included. As in Ref. |^0|], the solution is general 
in the sense that it expresses the full amplitude of the multiple processes of the Andreev 
electron ^ hole conversion and ordinary scattering via the amplitudes of the elementary 
processes. By this, the problem is split into independent and more simple problems. The 
theory of multi-channel bound states is also considered. The formulation which operates 
with excitations rather than bare particles, is especially convenient for the kinetic theory 
in the framework of the Boltzmann-type equations, for which it provides the boundary 



condition for the distribution function of the excitations [|Tl|,0. 

For a general case i. e. when the disorder and inelastic collisions are allowed, the boundary 
value of the 2-component wave functions = (") factorizing the trajectory Green's function 
are found in Sect.0. Since the mean field equations are linear, this result can be recast as 
the boundary condition for the Andreev amplitudes u/v of the Riccati equation approach. 
In a most compact and symmetric form, the boundary condition is given by Eq. (|5.5|) . For 
the specular interface, the boundary condition for the Riccati equation is given by Eq.(|5.7|), 
or explicitly by Eqs.( ^.8|) , and ( ^.9|) . 

One more form of the boundary condition is presented in Sect. |V B] , where the expression 
for the transfer matrix is derived. The transfer matrix, which couples the wave functions or 
the 1-point Green's functions on the chosen pair of in- and out-channels, absorbs information 
about all other 2(A^ — 1) channels. This modification of the boundary condition is convenient 
when one solves the Eilenberger equations for the 1-point Green's function. In the simplest 
2-channel case (specular reflection), this boundary condition reproduces Zaitsev's results 0. 
The new form seems to be more flexible and convenient. 

For the derivation, we use the technique of the 2-point Green's function. In our opinion, 
the technique provides an adequate language to discuss the semiclassical physics in super- 
conductors which we qualitatively considered in Section |. The 2-point Green's function 
gives a full description of the coherent propagation of electron and hole along a common 
classical path. In spite of the fact that observables can be expressed via the 1-point Green's 
function only, the language of the quasiclassical 2-point Green's function on is not redun- 
dant: Offering a physically transparent formalism, it is free from some uniqueness problems 
which plague the standard "^-integrated" formulation. Note also that with all possible sim- 
plifications already done, the quasiclassical 2-point Green's function obeys Eqs. (|2.2| ), and 
( ^.31 ) which, unlike the Eilenberger equations, have a familiar form of an equation for a 
propagator. Therefore, one may directly apply the intuition and experience gained in other 
fields of the quantum theory. 

Another attractive feature of the 2-point Green's function technique is that it allows 
one to define effective wave functions. The latter factorize the Green's function averaged 



32 



with respect to disorder or phonons. Although these "wave functions" have usual quantum 
mechanical meaning only in ballistic case, it seems to be advantageous that one may use the 
unified language of trajectories and wave functions discussing both the ballistic motion and 
the propagation in the presence of disorder or inelastic collisions. 

The effective wave function, = (") obeys the linear Andreev-type equation Eq.( |2.6| ). 
There is a variety of methods one can chose to solve the system of two linear differential 
equations for u and v. One of them is to derive the equation for the ratio u/v which turns 



out to be the Riccati equation suggested in p^ , p5| . As the logarithmic derivative ^p'/■^jJ 



in the usual Schrodinger equation does, the choice of the ratio u/v has the indisputable 
practical advantage which is due to insensitivity of the ratio to the normalisation of (p. The 
Riccati equation approach which has proven to be very convenient and efficient for numerics, 
finds rather natural physical interpretation in the 2-point Green's function technique of the 
present paper. (For the latest development of the Riccati equation approach including the 
interface boundary condition see e-preprint of M. Eschrig [P7|] .) 

An important part of this paper is the understanding that the classical trajectory trans- 
forms to a topologically 1-dimensional simply connected tree in the case of many interfaces 
and/or boundaries. The extended arguments in favour of this point of view have been pre- 
sented in Section |. Although, this assertion may look wrong in simple idealized geometries, 
like e.g., a sandwich with strictly parallel outer and the interface planes FigJ^, we argue 
that small deviations from the perfection eliminate accidental crossings of trajectories (as 
in non-integrable billiards). In our opinion, the difficulties with the quasiclassical theory 
encountered in |]^,|| are due to the fact that some interference contributions survive the 
procedure of the integration with respect to the layer thickness: Indeed, rigid variations of 
the layer thickness do not eliminate all the loops. We believe that some roughness, larger 
then the Fermi wave length but small and invisible on the quasiclassical scale, will restore 
the quasiclassical results. 

To show the new theory in action, we solve in Sect.|VI| two simple problems: (i) a film 
separated from a bulk material by a partially transparent interface; (ii) a two layer system 
with arbitrary transparent interface. (The latter was classified in [§] as quasiclassically 
unsolvable.) Motivated by recent ideas about the origin of the paramagnetic effect ||36|| , 
we evaluate the density of states and the superfluid density when the phase of the order 
parameters in the layers differs in vr, a scenario of paramagnetic instability suggested in 
Our results confirm the very possibility that the superfiuid density ps may be negative 



(Meissner "anti-screening") but we observe also that ps is strongly affected by reflection 
on the interface: when the probability of the reflection R > 0.04, the Meissner screening 
is restored. The implications of these results for a realistic theory of the paramagnetic 
instability requires further studies. 
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APPENDIX A: ADVANCED GREEN'S FUNCTIONS 



The advanced Green's function g^{xi,X2) is constructed in the same manner as the 
retarded one: One finds 0i,2 from Eg. (|2^) with substitutes for and builds the 
Green's function as in Eqs.(pT2D, ( ^J3D and ( ^ ). 

Combining Eqs.( |2.5| ) and ( p.7| ), one can see that the f^. [(p^) with cj)^ from Eg. ( p.6D 



satisfies the corresponding equation in the A-case. Then, the normahzed Eq.( 2.11 ) solutions 
are related to each other as 



^i = ^f44>fY , <pi = zf44>^y . (Ai) 

The Andreev amplitudes a and b Eq. (|2.16|) are related now as 



a 



^-l/(a^)* , 6^ = l/(6«)* 



and Green's functions as 

g^{xi,X2) = T^g'^* {xi,X2)% 

For future references, the symmetry in the 1-point Green's functions is given by the 
following well-known relations { £ = {e,n),£* = {e*, n)) 

9^{r) = -r. , g^{r) = {g^Ar))^ 



The first of them follows from Eq.(|Al|), and the second one refiects the symmetry = 

APPENDIX B: EVOLUTION IN PERIODIC POTENTIAL 

To prove validity of Eq . (12.261 ), one first solves the 2x2 eigenvalue problem 

Ul{x)iP{x) = 'yip{x) . 

and finds the eigenfunctions tpi^2 (with a; as a parameter) and the eigenvalues 71^2- It follows 
from the conservation of normalization in Eq. (|2.9| ) that Det U = 1, and, therefore, 

7i72 = 1 • 

Denote 71 the eigenvalue for which |7i| < 1, Q and normalize the eigenfunctions to satisfy 
ijj2'^i = 1- It is clear now, that ipiix) continued along the trajectory with the help of the 
evolution matrix Ul{x) gives the solution denoted in Eq. (|2.1CI| ) as (f>+(x): Indeed, it satisfies 
Eq.(p.6|) and decays at x — 00 as 'Ji^^ ■ By the same argument, 0_ = ip2- From Eq. ( p.20| ), 
the Green function now reads 



^ When considering g^, the variable e has a finite imaginary part and the matrix is not 
Hermitian. Then, the evolution matrix is not unitary and |7i,2| 7^ 1- 
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as 



Seeing that the evolution 2x2 matrix can be expanded in its normahzed eigenfunctions 

Ul{x) = ^(71 + 72) i + ^(71 - 72) + V'2V'l) , 

the traceless part of Ul{x) is proportional to g^. The normalization condition fixes the 
proportionality coefficient, and one comes to Eq.( p.26 ). 

To build the evolution matrix, one may use the following procedure. First consider 
two fundamental solutions to Eq. (^^) , ipj and ipii, which satisfy the following boundary 
conditions 

i^iix) = ( ) ' ^^^^^ + L) = J 
and find ipi{x + L) and ipii{x), 



+ L)= e^*-(^) ^^^^^ J , i^nix) = e^*-(^) ( 

The exponential factor is same for ipj and ipn as required by the conservation of the normal- 
ization in Eq.( |2.9D . The parameters a/, (a;), and ^l{x) can be calculated conveniently 
in the Riccati equation technique. 

Denote ao(a;; Xq) solution to Eq. (|2.23|) with the boundary condition aQ{x = Xq] Xq) = 0; 
then (see Eq.(|;2§) 



x+L 

aiix) = ao(x + L;x) , = 



Similarly, 

Piix) = po{x,x + L) 

where (/5o(x, xo)~ is the solution to Eq. (p.23|) with the boundary condition /5o( X — Xq, Xq J 
0. 

Building the evolution matrix from the fundamental solutions, one gets 



1 -I3dx) 



The traceless part of it, f/^(a;), reads 
& 

' 2 \2aL[x) -i 

Up to the normalization factor, this matrix is equal to g^{x 



' 2 \2aL{x) -1 + e-2^*^("') -ai(a;)/3i(a;) 
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APPENDIX C: FORMAL SOLUTION: I 



Here, we analyze some formal linear algebra aspects of the matching conditions in 
Generally, the wave functions may be presented in the following form: 

= ( ) ' ^^ = B^(^['^ (CI) 

Denote \X) the column with elements Xi, . . . ,Xn or Xii, . . . ,Xni. One obtains from 
Eqs.(|0|), and (g^) 

\A) = Su\B) , \B) = S^fi\A) . (C2) 

where fi and z> are diagonal matrices N x N with {fi)k'k' = l^k' and {P)kk = ^k- 
The two equalities in Eg . ( |C^) are compatible only if 

'^(M^ = det 1 - SuS^jl = , (C3) 



As expected, the parameters /I's and i/'s are not independent: Eq. (|C^) gives a relation 
among them, which is linear in each of the parameters (see Section |^), making it possible 
to express one of the /x's or z/'s through all others. 

For instance, one may give any values to /x's and z/'s in all channels excepting the /— th 
incoming one. Then ui is fixed by Eq.(^3|), 

= {l\S^i,Si\l) . (C4) 

where 

Si = (s^ -u^^^S^p)'^ , (C5) 



and z/*-'-* denotes the matrix which differs from z> only in that the element (z/*^'-*);/ = 0; here 
and below {i\Q\j) = Qij ■ 

From Eq . ( |C2|) one finds now the coefficients A's and B's, they are proportional to one of 
them, say Bi. It is convenient to put Bi = Cuf^. In this way, one gets the solution to the 
matching conditions in Eq. corresponding to the given set of fik' and z/j^;'s : 



here ui is given by Eq. (P^,A:, m = 1, . . . , N, and C is arbitrary. 

In the same way one builds the solution where all the /i's and z/'s are given as input 
excepting adjusted to meet the condition in Eq.( 



k',m = l,...,N, 

f,-^ = {n'\SuSl\n') . (C8) 



and 



sl=(S-fi^^'^SuY' , (C9) 
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APPENDIX D: FORMAL SOLUTION: II 

Transforming 1 - SuS^fi = i} ^ Su^^^S^ft^ (^1 - z/W (^S^ - u^^^S^fiY^ and 
using identity det(l + \X){Y\) = 1 + {Y\X), one gets 

Vi{u}, {^^}) = M) (l - z.«(/|5t^5H0) (Dl) 

with Si in Eg. (|C5|) , and 0^ denotes the matrix which differs from P only in that the element 
= 0. 
Similarly, 

Vi{u}, {^^}) = Vi{u}, (l - f^(-'Hn'\SuSl\n')) 

where S^^, is defined in Eq.(p9D. 

Sometimes calculations become shorter when one changes the representation. Since the 
5*- matrix in the Eq.( p.l|) is a scalar in the electron-hole space, the matching condition is 
unchanged by any rotation — ^ ipo = Oip, 

6= ' 

After the rotation, the basis wave function in Eg . (|Cl| ) has the same form with (...) {■■■)o 





fJ'k' — fJ'O 








1 - f^k'^0 












1 - /ioZ/i 


^Ok' 


1 - /ifc/Z/Q 

= Ay 




1 - /ioZ/Q 


Boi 


1 - /ioZ/i 
= ID, I 




1 - /ioZ/Q 



One sees that by proper rotations any pair /i^/, Vi can be nullified in the intermediate calcu- 
lations. Of course, all other coefficients will also be changed. Calculations done, one gets to 
the original basis. 

Egs.( [4.8|) , and (|4.9|) can be written in a more compact form. From Eg. (pl|) , 

1 ^ 'Di_ 
1 - at^af^ V, 

where 

Po = ^^({«^^)},{«^')}) 

and 

Pz = P({«(^)},{(a('^))«}) 
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Absorbing Vi into Si i.e. Si = ViSi, and using obvious = etc., etc the 

scattering amplitudes read 



Bi 



(0 



L({l\S^a^^)Si\l)-Via'i'^ 



4' = ^{k'\P^'^SiP^'^*\l) 



APPENDIX E: TRANSFER MATRIX 

Another possibihty for resolving the matching conditions is via the transfer matrix Ain'^i 

tpn' = Mn'^l^l , 

which couples the wave functions 

, f Ul\ . [ Un' 

\Vl J \ Vn' 

on a selected pair of trajectories / and n'; the parameters ^ik'^n' and Vi^i are supposed to be 
given. 

As usual, the transfer matrix can be built out of the elements of two particular solutions 
Take to be the solution in Eg. (|(J6D with yU„/ put to zero, 



and the solution Eq.(|C3) with i// = 0, 



1 

{l\S^fii^')Snn\l) 



where 

Sn'l 



Requiring that the transfer matrix reproduces the relations between ipn' and ipi in the two 
solutions, one gets the following result 

Mn'^i = {l\sUnr'{tc 1) (El) 

where 

A = {n'\S^^n\l){l\sUn') - {l\S^ fi'-^'''^ Sln\l){n'\Su^') S^n') , 
B = {n'\Su(^)sti\n') , C = {l\S^ fi^'^') Sln\l) ■ 
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The determinant of the transfer matrix is 



Det 



The inverse matrix reads 

M-l^ = {n'\S^n\l)-' 



{n'\Sn'i\i) 



1 -B 
C A 



(E2) 



Applying the matching conditions for the conjugated waves ip in Eq.( |3.3| ), one can check 
that the the corresponding transfer matrix is given by the inverse of that for ip i.e. 

Since ipn'^ipn'^ — i^l.^^ = '^z^VP'*) the conservation law in Eq.( p.9|) is not 
affected by knots. 



a. Transfer matrix 2x2 case 

The transfer matrix for the case when the knot mixes 2 in- to 2 out-trajectories, can be 
obtained from the general expression in Eq.(|El|). Algebraic simplifications of rather awkward 
expression gives a pretty compact result. Here, an alternative derivation, algebraically more 
transparent, is presented. 

Call the trajectories of interest by 1 and 1', and consider calculation of Air^i for given 
fi2' and 1^2 in Eq.(|Cl|). First note superconductivity influences the transfer matrix only via 
the trajectories 2 and 2'. Note also that in the normal metal case when /i2' = = 0, the 
transfer matrix is simply 

where s = Si>i. 

Since the 5— matrix is a electron- hole scalar, Eq.( |3.1| ) is invariant relative to rotations 
in the electron-hole space, ip ~^ "^o = C>ip, and one can resolve the matching conditions in 
arbitrary basis. 

The rotation 

a = - -''A ,d-' = (^ 

transforms (^^ ) to (^^,)^ = (J) and {^^) to (^^) as if in the normal state. Therefore, after 
the rotation, the transfer matrix is given by Eq. (p3|) , whereas in the original picture 

Inserting Eq.(^), the transfer matrix reads 
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One recognizes the combinations entering in Eq.( |2.17l ) with the important difference that 
i>2 and /i2' are parameters of the wave functions not at the same point but across the knot. 
Finally, 

-Mi'^i = ^ fl - T^92'.2] (E5) 



2s* V 1 + ^ 

where S = R = 1 — S, and 

_ 1 f ^ + f^2'i^2 -2z/2 

92'»2 



1 — /i2'Z/2 V 2/^2' — 1— /i2'Z^2 

is normalized ^|/,2 = 1, "across-the-knot" Green function. It can be presented in a factorized 
form as follows 



^ (1 + g2'.2) = - — - — ( M ( 1 , -^^2 ) 



(E6) 



I 1 \ (2') (2') 

Up to normalization, =0+ where 0Y is the knot value of 0+ on the trajectory 



/i2' 

2', and ( 1 , — z/2 ) = 0^'*- Taking into consideration Eqs.( |2.1^ ), and ( p.21| ), one concludes 
from Eq.(|E|) that (1 + ^2'.2) oc (1 + ^2') (1 + 92) so that 

(1 + ^2'.2) = AT-i (1 + ^2') (1+^2) , ^^= ^Sp (1 + 92') (1 + 92) , (E7) 

where the normalization |Sp (1 + ^2'»2) = 1 fixes the proportionality coefficient J\f. This 
formula expresses the "across-the-knot" function via Green's function on the trajectories 2 
and 2'. 

After some algebra, one gets another form of Eq.( 



fl'2'.2 = — TT^^ — r-r (92' +92 + 7;[92',92]- 



1 + 2192', 92]+ V 2 

(The anticommutator ^[§2', 92]+ in the denominator is proportional to the unit matrix and 
does not pose any problem). 
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FIGURES 




FIG. 1. Scattering on a partially transparent specular interface. The interface is depicted as the 
shaded region. The arrows show the direction of the (electron) velocity. The in-coming (out-going) 
trajectories are denoted 1 and 2 (1' and 2'). The filled circle, the knot (see text), is the "black 
box" where the scattering occurs. 
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FIG. 2. The typical trajectory in an ideal sandwich with the layers of an equal thickness and 
parallel surfaces. 
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FIG. 3. The typical trajectory in a sandwich with a rough surface. The roughness is shown 
schematicahy as a step, W being the height of the step. Unhke the ideal case in Fig. I, the paths 
1' and 2' return to the interface at different points. 




FIG. 4. An example of a tree-like trajectory. Pieces of the straight lines show the trajectories 
before or after they enter a knot (filled circles), i.e. before or after a collision with an inter- 
face. There is only one path connecting any two points on the tree so that the tree is effectively 
1-dimensional. 
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1 2' 

FIG. 5. Simplest 2x2 knot with two incoming 1 and 2, and two outgoing channels 1' and 2' 
(schematically) . 
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FIG. 6. For a many channel knot, one chooses a pair of trajectories, one in- and one out- 
(denoting them 1 and 1'), and considers them as a single trajectory with a knot on it. The transfer 
matrix relates to each other the wave functions across the knot. 
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(a) (b) 
FIG. 7. The typical trajectory formed by the total reflection on the outer surface and the partial 
reflection/transmissions on the interface (a). The numbers serve as markers for both direction and 
position. In (b), the structure of the tree-like trajectory is shown with the numbering as in (a). 
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(d) 

FIG. 8. The density of states (trajectory resolved) versus energy e/A at the interface of a bulk 
superconductor with the pair potential A; and a film A,. = —A; (see Fig. 0). The film thickness 
d is measured along the trajectory in units t'/||A|. The interface reflectivity R = 0, 0.1, 0.5, 0.9 in 
(a,b,c,d) respectively. 
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(a) (b) 

FIG. 9. (a )Real space classical trajectories of a particle in a two layers system formed by 
multiple reflections on the outer surface and the interface between layers . Numbers tag both the 
position of the particle on the trajectory as well as the direction, (b) The structure of the tree-like 
trajectory is shown. The points in real space and on the tree are marked by the same numbers in 
(a) and (b). 
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(a) 




(d) 

FIG. 10. Trajectory resolved density of states versus energy e/A at the interface. The order 
parameter A; = — A^. The thickness of the both layers is v/\/S.i\. The reflection i?, shown at the 
top of the pictures, is 0, 0.1, 0.6, and 0.9 in (a), (b), (c), and (d), respectively. 
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FIG. 11. The effective superfiuid density of a system of two layers, I and r, of equal thickness 
di = dr = v/\/S.i\ with the tt phase difference A/ = — for different reflectivity of the interface 
i? = 0, 0.01, 0.04 
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